use strict;
use warnings;
use Data::Dumper;
use Math::Combinatorics;

my @peptidos_formateados;
my %archivos_diferentes;
my (%hash_masas,%hash_intensidades,%hash_rt);
my $ruta_raw_file="G:\\00_0lymphos\\outs_antiguos\\IMAC_SCX1_081216";
my $fasta='G:\00_0lymphos\20080404uniprot_sprot_trembl_human.fasta';
my $path_DB="<$fasta";

my %nom_corto_largo;
my %proteinas_usadas;

my @mensajes_error;

open LISTA, "<IMAC_SCX1_081216.generar.out";
open OUT,">archivo.out";
open SALIDA,">resultados_comparados_Ascore.txt";
#open SAL, ">debug.txt";

principal();



sub principal{
	my $precision=0.5;
	my @datos_entrada;
	my @peptidos;
	my @pept_limpios;
	my @info_RAW;
	my @masas;
	my @intensidades;
	my @cargas;
	my @masas_precursor;
	my @scans;
	my @archivos;
	while(<LISTA>){
		chomp;
		push @datos_entrada,$_;
		my @d=split /\t/,$_;
		push @archivos,$d[0];
		push @scans,$d[1];
		my $pep=$d[7];
		$pep=(split /\./,$pep)[1];
		push @pept_limpios,$pep;
		$pep=~s/S\#/s/g;
		$pep=~s/T\#/t/g;
		$pep=~s/Y\#/y/g;
		$pep=~s/S\@/u/g;
		$pep=~s/T\@/v/g;
		$pep=~s/M\*/m/g;
		push @peptidos,$pep;
		push @cargas,$d[6];
		push @masas_precursor,($d[4]+$d[5]); 
	}
	for (my $i=0;$i<$#archivos+1;$i++){
		push @{$archivos_diferentes{$archivos[$i]}},$scans[$i];
	}
	#Obtengo masas, intensidades
		foreach my $raw(keys %archivos_diferentes){
			obtener_info_espectros($raw);
		}
	
	#Busqueda de peptidos en DB fasta
		foreach(@pept_limpios){
			$_=~s/\#//g;
			$_=~s/\@//g;
			$_=~s/\*//g;
		}
		
		my @resultados_find_proteins=find_proteins(\@pept_limpios,$path_DB);
			#RETURN: an array in the form of @array_returned=[@array of AC_nums matched,@array of protein sequences,
																#@array of peptide position inside protein].
		my @resultados_find_proteins_string;
		foreach my $res(@resultados_find_proteins){
			my @arr1= @{$res};
			my @arr2=@{$arr1[0]};
			my $txt="";
			foreach my $t(@arr2){
				$txt=$txt."||".$t;
				}
			push @resultados_find_proteins_string,$txt;	
			}
		print "Proteinas en memoria\n";
		
	#Calculo Ascores y genero outputs
	for (my $i=0;$i<$#peptidos+1;$i++){ 
		
		print "Calculado ascore de: ",$peptidos[$i],"\n";
		my $cadena_masas=$hash_masas{$archivos[$i]}{$scans[$i]};
		my $cadena_intensidades=$hash_intensidades{$archivos[$i]}{$scans[$i]};
		
		my @mas=split/\|/,$cadena_masas;
		my @int=split/\|/,$cadena_intensidades;
		shift @mas;
		shift @int;
		
		my @resultado=ascore($peptidos[$i],$precision,\@mas,\@int,$cargas[$i]);
		my $ascore_out;
		if($#resultado>7){#Peptidos con 2 fosforilaciones
			print SALIDA $datos_entrada[$i],"\t",$resultado[0],"\t",$resultado[1],"\t",$resultado[2],"\t",$resultado[3],"\t",$resultado[4],"\t",$resultado[5],"\t",$resultado[6],"\t",$resultado[7],"\t",$resultado[8],"\n";
			$ascore_out=$resultado[5]."//".$resultado[7];
		}
		elsif($#resultado==3){#Peptidos con el mismo numero de fosforilaciones que sitios fosforilables
			print SALIDA $datos_entrada[$i],"\t",$resultado[0],"\t",$resultado[1],"\t",$resultado[2],"\t\t\t",$resultado[3],"\n";
			$ascore_out=$resultado[3];
		}
		
		else{#Peptidos con 1 fosforilacion
			print SALIDA $datos_entrada[$i],"\t",$resultado[0],"\t",$resultado[1],"\t",$resultado[2],"\t",$resultado[3],"\t",$resultado[4],"\t",$resultado[5],"\t",$resultado[6],"\n";
			$ascore_out=$resultado[5];
		}
		
		my $ascores_reasignacion=0;#Testo si el peptido cambia de punto de fosforilacion tras calculo de Ascore
		unless($peptidos[$i] eq $resultado[0]){
			$ascores_reasignacion=1;
		}
		print OUT $datos_entrada[$i]."\t".
		$resultados_find_proteins_string[$i]."\t".
		$hash_masas{$archivos[$i]}{$scans[$i]}."\t".
		$hash_intensidades{$archivos[$i]}{$scans[$i]}."\t".
		$hash_rt{$archivos[$i]}{$scans[$i]}."\t".
		$ascore_out."\t".$resultado[0]."\t".$ascores_reasignacion."\n";
	}
	print OUT '###secuencias'."\n";
	foreach(keys %proteinas_usadas){
		print OUT $_."\t".$nom_corto_largo{$_}."\t".$proteinas_usadas{$_}."\n";
	}
}	
	
	
sub principal_dtas{
	my $precision=0.5;
	my @datos_entrada;
	my @peptidos;
	my @dtas;
	my @masas;
	my @intensidades;
	my @cargas;
	my @masas_precursor;
	while(<LISTA>){
		chomp;
		push @datos_entrada,$_;
		my @d=split /\t/,$_;
		my $pep=$d[2];
		$pep=~s/S\*/s/g;
		$pep=~s/T\*/t/g;
		$pep=~s/Y\*/y/g;
		$pep=~s/S\@/u/g;
		$pep=~s/T\@/v/g;
		push @peptidos,$pep;
		push @dtas,$d[6];
		my @datos_dta=lee_dta($d[6]);
		push @masas, $datos_dta[0];
		push @intensidades, $datos_dta[1];
		push @cargas,$datos_dta[2];
		push @masas_precursor,$datos_dta[3]; 
	}
	###
	
	
	for(my $i=0;$i<$#peptidos+1;$i++){
		print 1+$i,"\n";
		my @resultado=ascore($peptidos[$i],$precision,$masas[$i],$intensidades[$i],$cargas[$i]);
		
		if($#resultado>5){
			print SALIDA $datos_entrada[$i],"\t",$resultado[0],"\t",$resultado[1],"\t",$resultado[2],"\t",$resultado[3],"\t",$resultado[4],"\t",$resultado[5],"\t",$resultado[6]."\n";
		}
		else{
			print SALIDA $datos_entrada[$i],"\t",$resultado[0],"\t",$resultado[1],"\t",$resultado[2],"\t",$resultado[3],"\t",$resultado[4],"\n";
		}

	}
}

sub ascore{
	my $peptido=$_[0];			#en formato [styuv]
	my $precision=$_[1];		#en unidades m/z
	my @masas=@{$_[2]};			#referencia a array de masas
	my @intensidades=@{$_[3]};	#referencia a array intensidades
	my $z=$_[4];				#carga
	
	#Pre-procesamiento del espectro
	#my @eliminacion_masas= eliminacion_masas_precursor(\@masas,\@intensidades,$z,$masa_precursor_dta,$precision);
	#@masas=@{$eliminacion_masas[0]};
	#@intensidades=@{$eliminacion_masas[1]};


	#El numero de sitios de fosforilacion ha de ser mayor que el de 
	#fosforilaciones para que el Ascore tenga sentido. Si el Ascore no tiene sentido, lo devuelvo como 1000
		my $numero_sitios_fosforilacion=0;
		my $numero_fosforilaciones=0;
		my $peptido_conteo_modificaciones=$peptido;
		while($peptido_conteo_modificaciones=~s/[styuv]/x/){
			$numero_sitios_fosforilacion++;
			$numero_fosforilaciones++;
		}
		while($peptido_conteo_modificaciones=~s/[STY]/x/){
			$numero_sitios_fosforilacion++;
		}
		if($numero_sitios_fosforilacion==$numero_fosforilaciones){
			return ($peptido,0,0,1000);#(peptido,nivel-max,pepscore,ascore1)
		}
	
	
	if($numero_fosforilaciones==0){
		return ($peptido,0,0,1000);#(peptido,nivel-max,pepscore,ascore1) Por si se cuela uno sin fosforilar
	}
	elsif($numero_fosforilaciones==1){
		my @resultado=ascore_1_fosforilacion($peptido,$precision,\@masas,\@intensidades,$z);
		return @resultado;
	}
	elsif($numero_fosforilaciones==2){
		my @resultado=ascore_2_fosforilaciones($peptido,$precision,\@masas,\@intensidades,$z);
	}
	else{
		print "el numero maximo de fosforilaciones es 2\n";
		return ($peptido,0,0,1000);#(peptido,nivel-max,pepscore,ascore1) Para mas de 2 fosforilaciones
	}
		
	
	
}




sub ascore_1_fosforilacion{
	my $peptido=$_[0];			#en formato [styuv]
	my $precision=$_[1];		#en unidades m/z
	my @masas=@{$_[2]};			#referencia a array de masas
	my @intensidades=@{$_[3]};	#referencia a array intensidades
	my $z=$_[4];				#carga
	
	
	#Caso especial: numero de deshidrataciones es igual al numero de posibles deshidrataciones
		my $peptido_conteo_modificaciones=$peptido;
		my $numero_DH=0;
		my $numero_sitios_deshidratacion=0;
		while($peptido_conteo_modificaciones=~s/[uv]/x/){
			$numero_DH++;
			$numero_sitios_deshidratacion++;
		}
		while($peptido_conteo_modificaciones=~s/[STuv]/x/){
			$numero_sitios_deshidratacion++;
		}
		if($numero_DH>0){
			if($numero_sitios_deshidratacion==$numero_DH){
				return ($peptido,0,0,1000);#(peptido,nivel-max,pepscore,ascore1)
			}
		}
	
	
	my @peptidos_posibles=@{cadenitas($peptido)}; 
	
	
	#Obtencion de los Peptide Scores
	my %peptide_scores;
	my %datos_matches_10_niveles;
	foreach my $pep(@peptidos_posibles){
		my %fragmentos_teoricos= %{prever_fragmentos($pep,$z)};
		$datos_matches_10_niveles{$pep}=matches_10niveles_series_yb_completas(\%fragmentos_teoricos,\@masas,\@intensidades,$precision,$pep);
	}
	my %weigths=(1=>0.5, 2=>0.75, 3=>1, 4=>1, 5=>1, 6=>1, 7=>0.75, 8=>0.5, 9=>0.25,	10=>0.25);
	my $suma_niveles=7;	#suma de 0.5+0.75+1+1+1+1+0.75+0.5+0.25+0.25
	foreach my $pep(keys %datos_matches_10_niveles){
		my $suma_datos_matches=0;
		my $peptide_score=0;
		for (my $i=1;$i<11;$i++){
			$peptide_score=$peptide_score+($weigths{$i}*$datos_matches_10_niveles{$pep}{$i}{score}/$suma_niveles);
			$datos_matches_10_niveles{$pep}{$i}{score_ponderado}=($weigths{$i}*$datos_matches_10_niveles{$pep}{$i}{score}/$suma_niveles);
		}
		$peptide_scores{$pep}=$peptide_score;
	}
	
	
	#Ordeno peptidos segun Peptide score
	my @peptidos_con_sumas_scores_ordenadas=sort {$peptide_scores{$b} <=> $peptide_scores{$a}} keys %peptide_scores;
	my $p1=$peptidos_con_sumas_scores_ordenadas[0];
	my $p2=$peptidos_con_sumas_scores_ordenadas[1];
	my $p3=$peptidos_con_sumas_scores_ordenadas[2];

	
	#testeo si, en el caso de haber un tercer candidato, su peptide score es tan bueno como el segundo;
	if (defined $p3){ 			
		if((defined $peptide_scores{$p3})and($peptide_scores{$p2}==$peptide_scores{$p3})){
			#Si p2 tiene el mismo Peptide Score que p3, debo elegir uno para usar en el ascore 
			#print "pep2 ($p2) y pep 3 ($p3) son iguales ($peptide_scores{$p2}), respecto a $p1\n";
			###Determino cual de los peptidos p2 y p3 tiene la fosforilacion mas proxima a p1
			my @p_trozos_peptido_1= split //,$p1;
			my @p_trozos_peptido_2= split //,$p2;
			my @p_trozos_peptido_3= split //,$p3;
			my ($punto1,$punto2,$punto3);
			for(my $i=0;$i<$#p_trozos_peptido_1+1;$i++){if($p_trozos_peptido_1[$i]=~/[styuv]/){$punto1=1+$i;}}
			for(my $i=0;$i<$#p_trozos_peptido_2+1;$i++){if($p_trozos_peptido_2[$i]=~/[styuv]/){$punto2=1+$i;}}
			for(my $i=0;$i<$#p_trozos_peptido_3+1;$i++){if($p_trozos_peptido_3[$i]=~/[styuv]/){$punto3=1+$i;}}
			my $p2_usar; #uso el peptido con menor diferencia a punto1
			$punto2=abs($punto1-$punto2);
			$punto3=abs($punto1-$punto3);
			if($punto2<$punto3){
				$p2_usar=$p2;
			}
			else{
				$p2_usar=$p3;
			}
			$p2=$p2_usar;
		}
	}
	
	my @peptidos_ascore;
	push @peptidos_ascore,$p1;
	push @peptidos_ascore,$p2;
	
	
	#Obtengo nivel maxima diferencia entre Peptide Scores
	my $nivel_maxima_diferencia=0;
	my $difss=0;
	for(my $i=1;$i<11;$i++){
		my $dif= $datos_matches_10_niveles{$p1}{$i}{score}-$datos_matches_10_niveles{$p2}{$i}{score};
		if($dif>$difss){
			$difss=$dif;
			$nivel_maxima_diferencia=$i;
		}
	}
	if($nivel_maxima_diferencia==0){
		$nivel_maxima_diferencia=5;#Esto ocurre cuando el pep1 y el pep2 son identicos; uso el 5 de forma arbitraria 
	}
	
	my @p_trozos_peptido_1= split //,$p1;
	my @p_trozos_peptido_2= split //,$p2;
	my ($punto1,$punto2);
	for(my $i=0;$i<$#p_trozos_peptido_1+1;$i++){
		if($p_trozos_peptido_1[$i]=~/[styuv]/){
			$punto1=1+$i;
		}
	}
	for(my $i=0;$i<$#p_trozos_peptido_2+1;$i++){
		if($p_trozos_peptido_2[$i]=~/[styuv]/){
			$punto2=1+$i;
		}
	}
	if($punto1>$punto2){
		($punto1,$punto2)=($punto2,$punto1);
	}
	
	my %fragmentos_p1=%{prever_fragmentos($p1,$z)};
	%fragmentos_p1=%{tratamiento_fragmentos($p1,\%fragmentos_p1,$punto1,$punto2)};
	my %fragmentos_p2=%{prever_fragmentos($p2,$z)};
	%fragmentos_p2=%{tratamiento_fragmentos($p2,\%fragmentos_p2,$punto1,$punto2)};
	
	
	
	
	my %datos_ascore;
	$datos_ascore{$p1}=matches_ascore(\%fragmentos_p1,\@masas,\@intensidades,$precision,$p1);
	$datos_ascore{$p2}=matches_ascore(\%fragmentos_p2,\@masas,\@intensidades,$precision,$p2);
			
	my $ascore;
	
	
	my $matches1=1+$#{$datos_ascore{$p1}{$nivel_maxima_diferencia}{tags_match_intervalo}};
	my $matches2=1+$#{$datos_ascore{$p2}{$nivel_maxima_diferencia}{tags_match_intervalo}};
	my @array1=keys(%fragmentos_p1);
	my $size=1+$#array1;
	my $proporcion=$matches1."//".$size." (".$matches2.")";
	
	if((not defined $datos_ascore{$p1}{$nivel_maxima_diferencia}{score})and(not defined $datos_ascore{$p2}{$nivel_maxima_diferencia}{score})){
		return ($p1,$nivel_maxima_diferencia,$peptide_scores{$p1},0,$proporcion);
	}
	if((not defined $datos_ascore{$p1}{$nivel_maxima_diferencia}{score})){
		return ($p1,$nivel_maxima_diferencia,$peptide_scores{$p1},0,$proporcion);
	}
	$ascore=$datos_ascore{$p1}{$nivel_maxima_diferencia}{score}-$datos_ascore{$p2}{$nivel_maxima_diferencia}{score};	

	#Para ascores muy pequenyos, por problemas de redondeo, es posible que salga un numero negativo (-1*10^-10, por ejemplo). 
		if($ascore<0){
			print "Ascore de peptido $p1 ($ascore)redondeado a 0";
			$ascore=0;
		}
	
	return ($p1,$nivel_maxima_diferencia,$peptide_scores{$p1},$p2."(".$punto1."-".$punto2.")"," ",$ascore,$proporcion);


}

sub ascore_2_fosforilaciones{
	my $peptido=$_[0];			#en formato [styuv]
	my $precision=$_[1];		#en unidades m/z
	my @masas=@{$_[2]};			#referencia a array de masas
	my @intensidades=@{$_[3]};	#referencia a array intensidades
	my $z=$_[4];				#carga

	
	my @peptidos_posibles=@{cadenitas($peptido)}; 

	
	
	#Obtencion de los Peptide Scores
	my %peptide_scores;
	my %datos_matches_10_niveles;
	foreach my $pep(@peptidos_posibles){
		my %fragmentos_teoricos= %{prever_fragmentos($pep,$z)};
		$datos_matches_10_niveles{$pep}=matches_10niveles_series_yb_completas(\%fragmentos_teoricos,\@masas,\@intensidades,$precision,$pep);
	}
	my %weigths=(1=>0.5, 2=>0.75, 3=>1, 4=>1, 5=>1, 6=>1, 7=>0.75, 8=>0.5, 9=>0.25,	10=>0.25);
	my $suma_niveles=7;	#suma de 0.5+0.75+1+1+1+1+0.75+0.5+0.25+0.25
	foreach my $pep(keys %datos_matches_10_niveles){
		my $suma_datos_matches=0;
		my $peptide_score=0;
		for (my $i=1;$i<11;$i++){
			$peptide_score=$peptide_score+($weigths{$i}*$datos_matches_10_niveles{$pep}{$i}{score}/$suma_niveles);
			$datos_matches_10_niveles{$pep}{$i}{score_ponderado}=($weigths{$i}*$datos_matches_10_niveles{$pep}{$i}{score}/$suma_niveles);
		}
		$peptide_scores{$pep}=$peptide_score;
	}
	
	#Ordeno peptidos segun Peptide score
	my @peptidos_con_peptides_scores_ordenados=sort {$peptide_scores{$b} <=> $peptide_scores{$a}} keys %peptide_scores;
	my $p1=$peptidos_con_peptides_scores_ordenados[0];
	#Detecto puntos de modificacion
	my @p_trozos_peptido_1= split //,$p1;
	my $f1=0;
	my $f2=0;
	for(my $i=0;$i<$#p_trozos_peptido_1+1;$i++){
		if($p_trozos_peptido_1[$i]=~/[styuv]/){
			my $punto=1+$i;
			if($f1){
				$f2=$punto;
			}
			else{
				$f1=$punto;
			}
		}
	}
	
	#Encuentro Peptidos-2 para cada modificacion
	my $peptido2_f1;
	my $peptido2_f2;
	for (my $i=1;$i<$#peptidos_con_peptides_scores_ordenados+1;$i++){
		my @p_trozos_peptido_x=split //,$peptidos_con_peptides_scores_ordenados[$i];
		if(($p_trozos_peptido_x[(-1+$f1)]=~/[A-Z]/)and($p_trozos_peptido_x[(-1+$f2)]=~/[a-z]/)){
			$peptido2_f1=$peptidos_con_peptides_scores_ordenados[$i];
			last;
		}
	}
	for (my $i=1;$i<$#peptidos_con_peptides_scores_ordenados+1;$i++){
		my @p_trozos_peptido_x=split //,$peptidos_con_peptides_scores_ordenados[$i];
		if(($p_trozos_peptido_x[(-1+$f2)]=~/[A-Z]/)and($p_trozos_peptido_x[(-1+$f1)]=~/[a-z]/)){
			$peptido2_f2=$peptidos_con_peptides_scores_ordenados[$i];
			last;
		}
	}
		
	
	#Obtengo niveles maxima diferencia entre Peptide Scores (P1 respecto a peptido2_f1 y peptido2_f2)
	my $nivel_maxima_diferencia1=0;
	my $difss1=0;
	my $nivel_maxima_diferencia2=0;
	my $difss2=0;
	{
		for(my $i=1;$i<11;$i++){
			my $dif= $datos_matches_10_niveles{$p1}{$i}{score}-$datos_matches_10_niveles{$peptido2_f1}{$i}{score};
			if($dif>$difss1){
				$difss1=$dif;
				$nivel_maxima_diferencia1=$i;
			}
		}
		if($nivel_maxima_diferencia1==0){
			$nivel_maxima_diferencia1=5;#Esto ocurre cuando el pep1 y el pep2 son identicos; uso el 5 de forma arbitraria 
		}
		for(my $i=1;$i<11;$i++){
			my $dif= $datos_matches_10_niveles{$p1}{$i}{score}-$datos_matches_10_niveles{$peptido2_f2}{$i}{score};
			if($dif>$difss2){
				$difss2=$dif;
				$nivel_maxima_diferencia2=$i;
			}
		}
		if($nivel_maxima_diferencia2==0){
			$nivel_maxima_diferencia2=5;#Esto ocurre cuando los scores del pep1 y el pep2 son identicos; uso el 5 de forma arbitraria 
		}
	}
	
	
	#Obtengo fragmentos dentro de los intervalos adecuados
		#1#Encuentro puntos complementarios a las f1 y f2 en cada peptido
		my @p_trozos_peptido_2f1= split //,$peptido2_f1;
		my @p_trozos_peptido_2f2= split //,$peptido2_f2;
		my $aa_complement_1;
		for(my $i=0;$i<$#p_trozos_peptido_2f1+1;$i++){
			if($p_trozos_peptido_2f1[$i]=~/[a-z]/){
				unless((1+$i)==$f2){
					$aa_complement_1=1+$i;
				}
			}
		}
		my $aa_complement_2;
		for(my $i=0;$i<$#p_trozos_peptido_2f2+1;$i++){
			if($p_trozos_peptido_2f2[$i]=~/[a-z]/){
				unless((1+$i)==$f1){
					$aa_complement_2=1+$i;
				}
			}
		}
		
		#2#Ordeno los puntos que marcan los intervalos adecuados
			my ($punto11,$punto12,$punto21,$punto22);
			if($f1>$aa_complement_1){
				$punto11=$aa_complement_1;
				$punto12=$f1;
			}
			elsif($f1<$aa_complement_1){
				$punto12=$aa_complement_1;
				$punto11=$f1;
			}
			else{
				print "errorrrrrrrrrrrrrrr!!!!\n";
			}
			
			if($f2>$aa_complement_2){
				$punto21=$aa_complement_2;
				$punto22=$f2;
			}
			elsif($f2<$aa_complement_2){
				$punto22=$aa_complement_2;
				$punto21=$f2;
			}
			else{
				print "errorrrrrrrrrrrrrrr!!!!\n";
			}
		
		#3#Calculo fragmentos
			#Fragmentos para fosforilacion 1
			my %fragmentos_p1_1=%{prever_fragmentos($p1,$z)};
			%fragmentos_p1_1=%{tratamiento_fragmentos($p1,\%fragmentos_p1_1,$punto11,$punto12)};
			my %fragmentos_p1_2=%{prever_fragmentos($peptido2_f1,$z)};
			%fragmentos_p1_2=%{tratamiento_fragmentos($peptido2_f1,\%fragmentos_p1_2,$punto11,$punto12)};
			#Fragmentos para fosforilacion 2
			my %fragmentos_p2_1=%{prever_fragmentos($p1,$z)};
			%fragmentos_p2_1=%{tratamiento_fragmentos($p1,\%fragmentos_p2_1,$punto21,$punto22)};
			my %fragmentos_p2_2=%{prever_fragmentos($peptido2_f2,$z)};
			%fragmentos_p2_2=%{tratamiento_fragmentos($peptido2_f2,\%fragmentos_p2_2,$punto21,$punto22)};
		
		
		#4#Obtengo datos-ascore: con matches de los fragmentos site-determining
		my $datos_ascore_p1_1=matches_ascore(\%fragmentos_p1_1,\@masas,\@intensidades,$precision,$p1);
		my $datos_ascore_p1_2=matches_ascore(\%fragmentos_p1_2,\@masas,\@intensidades,$precision,$peptido2_f1);
				
		my $datos_ascore_p2_1=matches_ascore(\%fragmentos_p2_1,\@masas,\@intensidades,$precision,$p1);
		my $datos_ascore_p2_2=matches_ascore(\%fragmentos_p2_2,\@masas,\@intensidades,$precision,$peptido2_f2);
		
		
		#print Dumper $datos_ascore_p1_1;
		
		
		#print $#{${$datos_ascore_p1_1}{$nivel_maxima_diferencia1}{tags_match_intervalo}};
		#exit;
		my $matches1_1=1+$#{${$datos_ascore_p1_1}{$nivel_maxima_diferencia1}{tags_match_intervalo}};
		my $matches1_2=1+$#{${$datos_ascore_p1_2}{$nivel_maxima_diferencia1}{tags_match_intervalo}};
		my @array1=keys(%fragmentos_p1_1);
		my $size1=1+$#array1;
		my $proporcion1=$matches1_1."//".$size1." (".$matches1_2.")";
		
		my $matches2_1=1+$#{${$datos_ascore_p2_1}{$nivel_maxima_diferencia2}{tags_match_intervalo}};
		my $matches2_2=1+$#{${$datos_ascore_p2_2}{$nivel_maxima_diferencia2}{tags_match_intervalo}};
		my @array2=keys(%fragmentos_p2_1);
		my $size2=1+$#array2;
		my $proporcion2=$matches2_1."//".$size2." (".$matches2_2.")";
		

		
		my $ascore1=${$datos_ascore_p1_1}{$nivel_maxima_diferencia1}{score}-${$datos_ascore_p1_2}{$nivel_maxima_diferencia1}{score};	
		my $ascore2=${$datos_ascore_p2_1}{$nivel_maxima_diferencia2}{score}-${$datos_ascore_p2_2}{$nivel_maxima_diferencia2}{score};	
		
		#Para ascores muy pequenyos, por problemas de redondeo, es posible que salga un numero negativo (-1*10^-10, por ejemplo). 
		if($ascore1<0){
			print "Ascore1 de peptido $p1 ($ascore1)redondeado a 0";
			$ascore1=0;
		}
		if($ascore2<0){
			print "Ascore2 de peptido $p1 ($ascore2)redondeado a 0";
			$ascore2=0;
		}

		
		return ($p1,$nivel_maxima_diferencia1." | ".$nivel_maxima_diferencia2,$peptide_scores{$p1},$peptido2_f1."(".$punto11."-".$punto12.")",$peptido2_f2."(".$punto21."-".$punto22.")",$ascore1,$proporcion1,$ascore2,$proporcion2);
		

}



sub lee_dta{
	my $espectro= shift;
	my $masa_precursor_dta;
	my $z;
	my (@masas,@intensidades);
	open SPEC,"<$espectro";
	while(<SPEC>){
		chomp;
		($masa_precursor_dta,$z)=split / /,$_;
		last;
	}
	while(<SPEC>){
		chomp;
		unless($_=~/\d/){
			last;
		}
		my ($m,$i)=split / /,$_;
		push @masas,$m;
		push @intensidades,$i;
	}
	return (\@masas,\@intensidades,$z,$masa_precursor_dta);
	
}

sub prever_fragmentos{
	# Module to predict a spectra using a sequence. It uses the 20 aa and allows phosphorilations on STY (+80:s,t,y),
	#water loss on ST (-18:u,v), and methionine oxidation (+18:m). It produces a hash containing the {label}->mass.
	#It calculates the y and b series for the charge state, 
	#ARGUMENTS: there are 3 args:
	#       - $_[0]-> sequence
	#       - $_[1]-> charge state (1,2 or 3)
	#RETURN: a hash with {label}->mass.
		my %masas_aa = (
	A=>71.03711,V=>99.06841,L=>113.08406,I=>113.08406,P=>97.05276,M	=>131.0405,
	F=>147.0684,W=>186.07932,G=>57.02146,S=>87.03203,T=>101.04768,C=>160.03069,
	Y=>163.06332,N=>114.04293,Q=>128.05858,D=>115.02694,E=>129.04259,H=>137.05891,
	K=>128.09496,R=>156.1011,s=>166.99836,t=>181.01401,y=>243.02963,u=>69.02143,
	v=>83.03708,m=>147.03542);
	#Non standard aa's -> s=S+79.96633, t=T+79.96633, y=Y+80, u=S-18.0106, v=T-18.0106, m=M+15.99492
	#la C esta CAM!, la considero PTM fija


	my @aminoacidos=split //,$_[0];
	#La carga maxima permitida es 1 para z=1 o 2, 
	my $carga_maxima_permitida_fragmentos;
	if($_[1]==1){
		$carga_maxima_permitida_fragmentos=1;
	}
	elsif($_[1]==2){
		$carga_maxima_permitida_fragmentos=1;
	}
	else{
		$carga_maxima_permitida_fragmentos=$_[1]-1;
	}


	my $acumulado=0;
	my (@serie_y_1,@serie_b_1,@serie_y_2,@serie_b_2,@serie_y_3,@serie_b_3);

	for (my $i=0;$i<$#aminoacidos+1;$i++){
		$serie_b_1[$i]=$masas_aa{$aminoacidos[$i]}+$acumulado+1;
		$acumulado=$acumulado+$masas_aa{$aminoacidos[$i]};
		}
	my $M=$acumulado+19;
	for(my $i=0;$i<$#aminoacidos;$i++){
		$serie_y_1[$i]=$M-$serie_b_1[($#aminoacidos-$i)]+1;
		}
	$serie_y_1[$#aminoacidos]=$M-$serie_b_1[0]+1;
	if($_[1]>1){
		for (my $i=0;$i<$#aminoacidos+1;$i++){
			$serie_y_2[$i]=($serie_y_1[$i]+1)/2;
			if (defined $serie_b_1[$i]){
				$serie_b_2[$i]=($serie_b_1[$i]+1)/2;
				}
			}
		if($_[1]>2){
			for (my $i=0;$i<$#aminoacidos+1;$i++){
				$serie_y_3[$i]=($serie_y_1[$i]+1)/3;
				$serie_b_3[$i]=($serie_b_1[$i]+1)/3;
			}		
		}	    	
}


	my %fragmentos;

	##Prueba de fragmentos especiales
	#my $M_perdida98=($M+$_[1]-98)/$_[1];
	#$fragmentos{"esp1"}=sprintf("%.3f",$M_perdida98);	
	#$fragmentos{"esp1b"}=sprintf("%.3f",$M_perdida98/$_[1]);	
	#my $M_perdidaagua=($M+$_[1]-18)/$_[1];
	#$fragmentos{"esp2"}=sprintf("%.3f",$M_perdidaagua);
	#$fragmentos{"esp2b"}=sprintf("%.3f",$M_perdidaagua/$_[1]);


	for (my $i=0;$i<$#aminoacidos+1;$i++){
		my $j=$i+1;
		if ($i<($#aminoacidos)){											# El b(n) bo existe
			my $num=$j;
			if($num<10){$num='0'.$num};
			$fragmentos{"b$num+1"}=sprintf("%.3f",$serie_b_1[$i]);
			if(($_[1]>1)and($carga_maxima_permitida_fragmentos>1)){													# Si carga>1
				$fragmentos{"b$num+2"}=sprintf("%.3f",$serie_b_2[$i]);
				if(($_[1]>2)and($carga_maxima_permitida_fragmentos>2)){												# Si carga>2
					$fragmentos{"b$num+3"}=sprintf("%.3f",$serie_b_3[$i]);
				}
			}	
		}
		if(($i>0)){	
			my $num=$i;
			if($num<10){$num='0'.$num};														# El y(0) bo existe
			$fragmentos{"y$num+1"}=sprintf("%.3f",$serie_y_1[$i]);
			 if(($_[1]>1)and($carga_maxima_permitida_fragmentos>1)){													# Si carga>1
				$fragmentos{"y$num+2"}=sprintf("%.3f",$serie_y_2[$i]);
				if(($_[1]>2)and($carga_maxima_permitida_fragmentos>2)){												# Si carga>2
					$fragmentos{"y$num+3"}=sprintf("%.3f",$serie_y_3[$i]);	
				}	
			 }
		}	
	}
		#Solo se consideran los fragmentos menos de m/z=2000
		foreach my $frag(keys %fragmentos){
			if($fragmentos{$frag}>=2000){
				delete $fragmentos{$frag};
			}
		}


	return \%fragmentos;
}

sub matches_10niveles_series_yb_completas{
	my %fragmentos=%{$_[0]};
	my @masas=@{$_[1]};
	my @intensidades=@{$_[2]};
	my $fragment_error=$_[3];
	my $peptido=$_[4];
	
	my $numero_fragmentos=0;
	foreach(keys %fragmentos){
		$numero_fragmentos++;
		#print $_,"\t",$fragmentos{$_},"\n";
	}
	
	
	my $num_picos_cada_100;
	my %datos_numeros_picos;
	my @rangos_masas;
	my ($rango_masas_inf,$rango_masas_sup)=(200,2000);
	for(my $i=$rango_masas_inf;$i<$rango_masas_sup;$i=$i+100){
		push @rangos_masas,$i;
	}
	for($num_picos_cada_100=1;$num_picos_cada_100<11;$num_picos_cada_100++){
		my @masas_exp_match_intervalo;
		my @masas_teor_match_intervalo;
		my @tags_match_intervalo;
		my $correspondencias=0;
		my $num_picos_dentro_intervalo=0;
		for my $intervalo(@rangos_masas){
			my %picos_intervalo;
			my $intervalo_inferior=$intervalo-100;
			for(my $i=0;$i<$#masas+1;$i++){
				if(($masas[$i]>$intervalo_inferior)and($masas[$i]<=$intervalo)){
					$picos_intervalo{$masas[$i]}=$intensidades[$i];
					$num_picos_dentro_intervalo++;
				}
			}
			my $num_tests_dentro_intervalo=0;
			my @masas_ordenadas_por_intensidad;
			
			if($num_picos_dentro_intervalo==0){#ningun pico experimental en el intervalo
										       #No hay picos en el intervalo\n";
				next;
			}
			elsif($num_picos_dentro_intervalo==1){ #un solo pico experimental en el intervalo
				foreach my $key(keys %picos_intervalo ){
					push @masas_ordenadas_por_intensidad,$key;
				}
			}
			else{#mas de 1 pico experimental en el intervalo
				foreach my $k (sort { $picos_intervalo{$b} <=> $picos_intervalo {$a}} keys %picos_intervalo ){
					push @masas_ordenadas_por_intensidad,$k;
				}
				
			}
			
			#	Aqui se realiza el matching entre masas teoricas y experimentales
			my $masas_testadas=0;
			foreach my $masas_exp(@masas_ordenadas_por_intensidad){
				$masas_testadas++;
				foreach(keys %fragmentos){
					if(abs($fragmentos{$_}-$masas_exp)<=$fragment_error){
						$correspondencias++;
						push @masas_exp_match_intervalo,$masas_exp;
						push @masas_teor_match_intervalo,$fragmentos{$_};
						push @tags_match_intervalo, $_;
					}
				}
				if($masas_testadas==$num_picos_cada_100){
					last;
				}
			}
		}
				#$numero_fragmentos=> trials
				#$correspondencias=>  successes
				#$prob=> numero picos cada 100 / 100

		my $probablity=($num_picos_cada_100/100);
		my $binomial_acum_rango=binomial($numero_fragmentos,$correspondencias,$probablity);
		my $score;
		if($binomial_acum_rango>0){
			$score=-10*(log($binomial_acum_rango)/log(10));
		}
		else{
			$score=0;
		}
		
		$datos_numeros_picos{$num_picos_cada_100}{masas_exp_match_intervalo}=\@masas_exp_match_intervalo;
		$datos_numeros_picos{$num_picos_cada_100}{masas_teor_match_intervalo}=\@masas_teor_match_intervalo;
		$datos_numeros_picos{$num_picos_cada_100}{tags_match_intervalo}=\@tags_match_intervalo;
		$datos_numeros_picos{$num_picos_cada_100}{score}=$score;
	}
	
	return \%datos_numeros_picos;
	


}

sub cadenitas{
	#This module uses only one argument: a peptide->$_[0]='AARTDRTsERuS'
	#It returns all possible combinations whith phosphorilation and
	#dehidratation sites constant (it uses combinations of 'n' possible sites taken in groups of 
	# 't' real sites without repetition). This combinatoric operation is used for phosphotylations
	#and dehidratations concurrently 
	my $cadena=$_[0];
	my @puntos_modificacion;
	my @puntos_deshidratacion;
	my @posibles_puntos_modificacion;
	my @aminoacidos= split //,$cadena;
	for(my $i=0;$i<$#aminoacidos+1;$i++){
		if(($aminoacidos[$i] eq 's')or($aminoacidos[$i] eq 't')or($aminoacidos[$i] eq 'y')){
			push @puntos_modificacion, $i;
			}
		if($aminoacidos[$i]=~m/[s,t,S,T,u,v,y,Y]/){
			push @posibles_puntos_modificacion, $i;
			}
		if($aminoacidos[$i]=~m/[u,v]/){
			push @puntos_deshidratacion, $i;
			}			
		}
	my $numero_modificaciones=@puntos_modificacion;
	my $posiciones_dh= @puntos_deshidratacion;		
	my $cadena_sin_modificaciones= $cadena;
	$cadena_sin_modificaciones=~s/s/S/g;
	$cadena_sin_modificaciones=~s/t/T/g;
	$cadena_sin_modificaciones=~s/u/S/g;
	$cadena_sin_modificaciones=~s/v/T/g;
	$cadena_sin_modificaciones=~s/y/Y/g;
		
	my %posibles_secuencias;
	my $combinat = Math::Combinatorics->new(count => $numero_modificaciones,data => [@posibles_puntos_modificacion],);
	while(my @combo = $combinat->next_combination){
		my @posiciones_modific_tras_fosfo;
		foreach my $elim (@combo){
			foreach(@posibles_puntos_modificacion){
				unless($elim==$_){
					push @posiciones_modific_tras_fosfo, $_;
					}
				}
			}
			my @pep_printar=split //,$cadena_sin_modificaciones;
			for(@combo){
				if($pep_printar[$_] eq 'S'){
					$pep_printar[$_]='s';
					}
				elsif($pep_printar[$_] eq 'T'){
					$pep_printar[$_]='t';
					}
				elsif($pep_printar[$_] eq 'Y'){
					$pep_printar[$_]='y';
					}	
				}
			my $peptido_printar= join('',@pep_printar);		
		if ($posiciones_dh>0){#si hay deshidrataciones 	
			my $combinat2= Math::Combinatorics->new(count => $posiciones_dh,data => [@posiciones_modific_tras_fosfo],);
			while(my @combo2 = $combinat2->next_combination){
				my @pep_printar2=split //,$peptido_printar;
				for(@combo2){
					if($pep_printar2[$_] eq 'S'){
						$pep_printar2[$_]='u';
						}
					elsif($pep_printar2[$_] eq 'T'){
						$pep_printar2[$_]='v';
						}	
				}
				my $peptido_printar2= join('',@pep_printar2);
				$posibles_secuencias{$peptido_printar2}=1; #Importante! Pongo el peptido en variable tras fosfo+deshidratacion
				}
			}
		else{#esto debo hacerlo solo si no hay deshidrataciones
			$posibles_secuencias{$peptido_printar}=1;#Importante! Pongo el peptido en variable tras 
													#fosfo solamente (no hay DH)	
		}
  }
  if ($numero_modificaciones==0){#caso especifico en el que solo haya deshidrataciones
  	if ($posiciones_dh>0){#si hay deshidrataciones
  		my $combinat2= Math::Combinatorics->new(count => $posiciones_dh,data => [@posibles_puntos_modificacion],);
  		while(my @combo2 = $combinat2->next_combination){
				my @pep_printar2=split //,$cadena_sin_modificaciones;
				for(@combo2){
					if($pep_printar2[$_] eq 'S'){
						$pep_printar2[$_]='u';
						}
					elsif($pep_printar2[$_] eq 'T'){
						$pep_printar2[$_]='v';
						}	
				}
				my $peptido_printar2= join('',@pep_printar2);
				$posibles_secuencias{$peptido_printar2}=1; #Importante! Pongo el peptido en variable tras fosfo+deshidratacion
		}
  	}
  }
	my @cadenitas;
	my $cadena_conteo_modificaciones=$cadena;
	my @cad=split //,$cadena;
	my $num_total_modif=0;
	foreach my $c(@cad){
		if($c=~/[styuv]/){
			$num_total_modif++;
		}
	}
	
	foreach my $seq(keys %posibles_secuencias){
		my @cad=split //,$seq;#esto lo hago para asegurarme que no se cuela ninguna secuencia erronea
		my $num_modif=0;
		foreach my $c(@cad){
			if($c=~/[styuv]/){
				$num_modif++;
			}
		}
		if($num_total_modif==$num_modif){
			push @cadenitas, $seq;
		}
	}
	return \@cadenitas;		
}

sub binomial{
	#Binomial distribution
	use Math::NumberCruncher;
	#($n,$k,$p)->$k successes in $n tries, given a probability of $p	
	my $n=$_[0];
	my $k=$_[1];
	my $p=$_[2];
	my $binomial_acumulada=0;
	for(my $i=$k; $i<$n+1;$i++){
		my $binomial=Math::NumberCruncher::Binomial($n,$i,$p);
		$binomial_acumulada= $binomial_acumulada+$binomial;
	}
	return $binomial_acumulada;
}

sub eliminacion_masas_precursor{
	my @masas= @{$_[0]};
	my @intensidades= @{$_[1]};
	my @masas_a_eliminar;
	my $z=$_[2];
	my $masa_precursor_dta=$_[3];
	my $precision=$_[4];
	for (my $i=$z;$i>0;$i--){
		push @masas_a_eliminar,($masa_precursor_dta+$i-98)/$i;#perdida 98
		push @masas_a_eliminar,($masa_precursor_dta+$i-18)/$i;#perdida agua
		push @masas_a_eliminar,($masa_precursor_dta+$i)/$i;#perdida 98
	}
	
	my @masas_procesadas;
	my @intensidades_procesadas;
	for(my $i=0;$i<$#masas+1;$i++){
		my $elimino=0;
		foreach my $me(@masas_a_eliminar){
			if(abs($masas[$i]-$me)<=$precision){
				print "Masa $masas[$i] ($intensidades[$i]) eliminada\n";
				$elimino=1;
			}
		}
		unless($elimino){
			push @masas_procesadas,$masas[$i];
			push @intensidades_procesadas,$intensidades[$i];
		}
	}
	return (\@masas_procesadas,\@intensidades_procesadas);

}

sub tratamiento_fragmentos{
	#argumentos-> $_[0]:%fragmentos,$_[1]:$peptido,$_[2]:primer aa intervalo, $_[3]: ultimo aa intervalo;
	########		Tratamiento fragmentos		######
	#Aqui genero un nuevo hash para poner los fragmentos que estan entre el primero y ultimo aa que se pueda fosforilar
		
		my $peptido=$_[0];
		my %fragmentos=%{$_[1]};
		my $primer_aa_fosforilable=$_[2];
		my $ultimo_aa_fosforilable=$_[3];
		my %fragmentos_filtrados;
	
		my $largo_pep= length $peptido;
		
		foreach my $frag (keys %fragmentos){
			my $frag2=(split /\+/,$frag)[0];
			if($frag2=~m/y/){
				$frag2=~s/y//;		#para las y, valen los fragmentos entre [largo-ultimo+1 , largo-primero]
				if(($frag2>=($largo_pep-$ultimo_aa_fosforilable+1))and($frag2<=($largo_pep-$primer_aa_fosforilable))){
					$fragmentos_filtrados{$frag}=$fragmentos{$frag};
				}
			}
			elsif($frag2=~m/b/){	#para las b, valen los fragmentos entre [primero , ultimo-1]
				$frag2=~s/b//;
				if(($frag2>=($primer_aa_fosforilable))and($frag2<=($ultimo_aa_fosforilable-1))){
					$fragmentos_filtrados{$frag}=$fragmentos{$frag};
				}
			}
		}
		
		
		%fragmentos=();
		%fragmentos = %fragmentos_filtrados;#Uso solo fragmentos en zona de interes		
	return \%fragmentos;
}


sub matches_ascore{
	my %fragmentos=%{$_[0]};
	my @masas=@{$_[1]};
	my @intensidades=@{$_[2]};
	my $fragment_error=$_[3];
	my $peptido=$_[4];
	
	
	my $numero_fragmentos=0;
	foreach(keys %fragmentos){
		$numero_fragmentos++;
	}
	
	my $num_picos_cada_100;
	my %datos_numeros_picos;
	my @rangos_masas;
	my ($rango_masas_inf,$rango_masas_sup)=(200,2000);
	for(my $i=$rango_masas_inf;$i<$rango_masas_sup;$i=$i+100){
		push @rangos_masas,$i;
	}
	for($num_picos_cada_100=1;$num_picos_cada_100<11;$num_picos_cada_100++){
		my @masas_exp_match_intervalo;
		my @masas_teor_match_intervalo;
		my @tags_match_intervalo;
		my $correspondencias=0;
		my $num_picos_dentro_intervalo=0;
		for my $intervalo(@rangos_masas){
			my %picos_intervalo;
			my $intervalo_inferior=$intervalo-100;
			for(my $i=0;$i<$#masas+1;$i++){
				if(($masas[$i]>$intervalo_inferior)and($masas[$i]<=$intervalo)){
					$picos_intervalo{$masas[$i]}=$intensidades[$i];
					$num_picos_dentro_intervalo++;
				}
			}
			my $num_tests_dentro_intervalo=0;
			my @masas_ordenadas_por_intensidad;
			
			if($num_picos_dentro_intervalo==0){#ningun pico experimental en el intervalo
										       #No hay picos en el intervalo\n";
				next;
			}
			elsif($num_picos_dentro_intervalo==1){ #un solo pico experimental en el intervalo
				foreach my $key(keys %picos_intervalo ){
					push @masas_ordenadas_por_intensidad,$key;
				}
			}
			else{#mas de 1 pico experimental en el intervalo
				foreach my $k (sort { $picos_intervalo{$b} <=> $picos_intervalo {$a}} keys %picos_intervalo ){
					push @masas_ordenadas_por_intensidad,$k;
				}
				
			}
			
			#	Aqui se realiza el matching entre masas teoricas y experimentales
			my $masas_testadas=0;
			foreach my $masas_exp(@masas_ordenadas_por_intensidad){
				$masas_testadas++;
				foreach(keys %fragmentos){
					if(abs($fragmentos{$_}-$masas_exp)<=$fragment_error){
						$correspondencias++;
						push @masas_exp_match_intervalo,$masas_exp;
						push @masas_teor_match_intervalo,$fragmentos{$_};
						push @tags_match_intervalo, $_;
					}
				}
				if($masas_testadas==$num_picos_cada_100){
					last;
				}
			}
		}

		my $probablity=($num_picos_cada_100/100);
		my $binomial_acum_rango=binomial($numero_fragmentos,$correspondencias,$probablity);
		my $score;
		if($binomial_acum_rango>0){
			$score=-10*(log($binomial_acum_rango)/log(10));
		}
		else{
			$score=0;
		}
		
		$datos_numeros_picos{$num_picos_cada_100}{masas_exp_match_intervalo}=\@masas_exp_match_intervalo;
		$datos_numeros_picos{$num_picos_cada_100}{masas_teor_match_intervalo}=\@masas_teor_match_intervalo;
		$datos_numeros_picos{$num_picos_cada_100}{tags_match_intervalo}=\@tags_match_intervalo;
		$datos_numeros_picos{$num_picos_cada_100}{score}=$score;
	}
	
	return \%datos_numeros_picos;



}


sub mayusc_minusc{
	#devuelve 1 si los 2 aa son minusculas y 0 si no
	my $letra1=$_[0];
	my $letra2=$_[1];
	
	if($letra1=~/[a-z]/){
		if($letra2=~/[a-z]/){
			return 1;
		}
	}
	return 0;
}


sub obtener_info_espectros{
	#obtengo espectros y Rt de cada scan de cada archivo RAW; los pongo en hashes del tipo $hash{$raw_file}{$scan}
	#		(%hash_masas,%hash_intensidades,%hash_rt)
	#	nota: aqui no importa que para un mismo scan_number tenga las dos alternativas de carga: tanto el espectro como 
	#		el tiempo de retencion son equivalentes
	my $raw= shift;
	my @array_scans=@{$archivos_diferentes{$raw}};
		my %resultado=obtener_espectro($ruta_raw_file."\\".$raw,\@array_scans);
		my %hash_rt_local= obtener_rt($ruta_raw_file."\\".$raw);
		foreach my $sc(keys %resultado){
			my @arrar_masas=@{${$resultado{$sc}}[0]};
			my @arrar_intensidades=@{${$resultado{$sc}}[1]};
			my $cadena_arrar_masas="";
			my $cadena_arrar_intensidades="";
			foreach(@arrar_masas){
				$cadena_arrar_masas=$cadena_arrar_masas."|".$_;
				}
			foreach(@arrar_intensidades){
				$cadena_arrar_intensidades =$cadena_arrar_intensidades."|".$_;
				}
			$hash_masas{$raw}{$sc}=$cadena_arrar_masas;	
			$hash_intensidades{$raw}{$sc}=$cadena_arrar_intensidades;
			$hash_rt{$raw}{$sc}=$hash_rt_local{$sc};	
		}
}

sub obtener_espectro{
    my $rawfile=$_[0];
    my @scans=@{$_[1]};
    push (@mensajes_error, $rawfile." no es un archivo\n") and return unless -f $rawfile; 
    use Win32::OLE;
    use Win32::OLE::Variant;
    my %resultados;
	Win32::OLE->Option( Warn => 3 );
	my $y = eval{
			Win32::OLE->new('XRaw.XRaw.1', sub{ $_[0]->Close } );
	};
	$y->Open( $rawfile );
	my $spectra= $y->Detector(0,1)->Spectra(0);
	
	foreach my $scan(@scans){
			eval{$spectra->Item($scan)->Data} or push @mensajes_error,"Espectro $scan falla\n";
			if($@){ 
				next;
			}
			my $espectro=$spectra->Item($scan)->Data;
			my @masas= @{@{$espectro}[0]};
			my @intensidades= @{@{$espectro}[1]};
			my %par_int_masa;
			for (my $i=0;$i<$#masas+1;$i++){
				$par_int_masa{sprintf("%.3f",$intensidades[$i])}=sprintf("%.3f",$masas[$i]);
				}
			my (@masas_ord,@int_ord);my $indice=0;
			foreach my $food (sort { $b <=> $a } keys %par_int_masa) {
				$masas_ord[$indice]=$par_int_masa{$food};
				$int_ord[$indice]=$food;
				$indice++;
				if($indice>1000){
					last; #IMPORTANTE: solo me quedo con las 1000 masas mas intensas
					}
			}        	
		   $resultados{$scan}=[\@masas_ord,\@int_ord];	
	}
	return %resultados;
}
   
sub obtener_rt{
	my $archivo_raw=$_[0];
	my %rt;
	use Win32::OLE;
	use Win32::OLE::Variant;

	my $enlaceRaw =Win32::OLE->new('XRawfile.XRawfile.1', sub{ $_[0]->Close } );
	$enlaceRaw->Open($archivo_raw);
	$enlaceRaw->SetCurrentController( 0, 1 );

	my $first = Win32::OLE::Variant->new(VT_I4|VT_BYREF,0);
	my $last = Win32::OLE::Variant->new(VT_I4|VT_BYREF,0);

		$enlaceRaw->GetFirstSpectrumNumber($first);
		$enlaceRaw->GetLastSpectrumNumber($last);
			
	my $cuenta_scans=1;
		for my $i ( $first .. $last ){
			# do not pass an initialization value to new or it won't work
			# this requires Activestate perl 5.8.*
			my $filter = Win32::OLE::Variant->new(VT_BSTR|VT_BYREF);
			my $t = Win32::OLE::Variant->new(VT_R8|VT_BYREF,0);
			$enlaceRaw->GetFilterForScanNum($i,$filter);
			$enlaceRaw->RTFromScanNum($i,$t);
			$t= $t*60;
			$rt{$cuenta_scans}=sprintf("%.3f",$t);        
			$cuenta_scans++;
			
		}
	return %rt;
}	
	
	
sub find_proteins {

# Module to extract find all proteins that match for a certain peptide inside a database. It returns
#an array with the values of the AC_Num (matching proteins) and sequence of the protein
#ARGUMENTS: there are 2 args:
#       - $_[0]-> reference to an array containing the peptides we want to evaluate (eg, find in what proteins they
                        # can be present)
#       - $_[1]-> Path to the database file in fasta format                        
#RETURN: an array in the form of @array_returned=[@array of AC_nums matched,@array of protein sequences,
                                                        #@array of peptide position inside protein].
#If no protein contains the searched peptide,
# an empty string is reurnd for this  $array_returned[$i] value
        
        my @peptides=@{$_[0]};
        my $database=$_[1];
                    
        my $temp="";
        my $nom="";
        my $nom_largo="";
        my @nombres;
        my @nombres_largos;
        my @secuencias;
        my %nom_seq;
        
        
        open PROTS,$database;

        while (<PROTS>){        #Entire fasta database is loaded into RAM. @nombres are the parsed
                if($_=~/^>/){   # ACNums, and @secuencias the sequences
                    $nom=$_;
                    $nom=~s/\>//;
                    $nom_largo=$nom;
                    ($nom,$nom_largo)=split / /,$nom,2;
                    chomp $nom;
                    chomp $nom_largo;
                    push @nombres,$nom;
                    push @nombres_largos,$nom_largo;
                    unless($temp eq ""){
                        push @secuencias,$temp;
                        }
                    $temp="";
                }
                else{
                    chomp;
                    $temp=$temp.$_;
                }
        }

        push @secuencias,$temp; #I must push the last sequence, that is still in $temp variable
        close PROTS;
        
        for(my $i=0;$i<$#nombres+1;$i++){       #A hash is built, using ACNum as keys, sequences as content
            $nom_seq{$nombres[$i]}=$secuencias[$i];     
        	$nom_corto_largo{$nombres[$i]}=$nombres_largos[$i];
        }
    	
    	
        my %pep_problema_prots;
        my @ac_nums;
        for(my $i=0;$i<$#peptides+1;$i++){      #The presence of the peptide sequence is searched into
            my @ac_nums_n;my @ac_nombres_largos; my @secs_n; my @posicion_primer_aa_n;          # all the protein sequences.
                for my $nomm(keys %nom_seq){
                    if((index $nom_seq{$nomm},$peptides[$i])>-1){# If a match is found, values of ACnum amb protein seq,
                        push @ac_nums_n,$nomm;                   # are pushed into arrays @ac_nums_n and @secs_n. The position
                        push @ac_nombres_largos,$nom_corto_largo{$nomm};#of the peptide inside protein is also found (first peptide's aa).
                        push @secs_n,$nom_seq{$nomm};           
                        push @posicion_primer_aa_n,(1+index $nom_seq{$nomm},$peptides[$i]);
                        $proteinas_usadas{$nomm}=$nom_seq{$nomm};
                    }
                }
                if ((defined $ac_nums_n[0]) and(defined $secs_n[0])){
                        $ac_nums[$i]=[\@ac_nums_n,\@ac_nombres_largos,\@secs_n,\@posicion_primer_aa_n];
                                #If arrays @ac_nums_n and @secs_n are defined,
                                #(at least one value->$ac_nums_n[0] and $secs_n[0])
                }               #they are referenced into the general array @ac_nums[$i] 
                else{
                        $ac_nums[$i]=[\@ac_nums_n,\@ac_nombres_largos,\@secs_n,\@posicion_primer_aa_n];
                        #If arrays @ac_nums_n and @secs_n are not defined, an emty string
                }       # is placed into the $ac_nums[$i] array
                		#
                
        }                                             

        return @ac_nums;
}		
	
	
	
